#pragma once
#include <vector>

using namespace std;
#ifndef Pi 
#define Pi 3.141592653589793238462643 
#endif 
 static class Binomial
{
public:
	static double generic_payoff(const double& S, const double& K,char PutCall, char OpStyle)
	{

		if(OpStyle == 'a')
		{
			if(PutCall == 'c')
				return max(0.0,S-K);
			else
				return max(0.0,K-S);
		}
		else
		{
			if(PutCall == 'c')
				return max(0.0,S-K);
			else
				return max(0.0,K-S);
		}

	}
	static double BlackScholes(char CallPutFlag, double S, double X, double T, double r, double v)
	{
		double d1, d2;


		d1=(log(S/X)+(r+v*v/2)*T)/(v*sqrt(T));
		d2=d1-v*sqrt(T);

		if(CallPutFlag == 'c')
			return S *CND(d1)-X * exp(-r*T)*CND(d2);
		else if(CallPutFlag == 'p')
			return X * exp(-r * T) * CND(-d2) - S * CND(-d1);
	}

	static double option_price_delta_american_binomial (char PutCall, char OpStyle, const double& S,const double& K, 
		  const double& r, 
		const double& sigma, const double& t, const int& no_steps)
	{
		double R = exp(r*(t/no_steps));
		double Rinv = 1.0/R;
		double u = exp(sigma*sqrt(t/no_steps));
		double d = 1.0/u;
		double uu= u*u;
		double pUp = (R-d)/(u-d);
		double pDown = 1.0 - pUp;
		vector<double> prices (no_steps+1);
		prices[0] = S*pow(d, no_steps);
		for (int i=1; i<=no_steps; ++i) prices[i] = uu*prices[i-1];
		vector<double> values (no_steps+1);
		for (int i=0; i<=no_steps; ++i) values[i] = generic_payoff(prices[i],K, PutCall,  OpStyle);
		for (int CurrStep=no_steps-1 ; CurrStep>=1; --CurrStep) {
			for (int i=0; i<=CurrStep; ++i) {
				prices[i] = d*prices[i+1];
				values[i] = (pDown*values[i]+pUp*values[i+1])*Rinv;
				values[i] = max(values[i], generic_payoff(prices[i],K, PutCall,  OpStyle));
			};
		}; 
		double delta = (values[1]-values[0])/(S*u-S*d); 
		return delta;
	};
	
	static double option_price_delta_european_binomial(char PutCall, char OpStyle, const double& S, const double& K, 
		 const double& r, 
		const double& sigma, const double& t, const int& no_steps)
	{
		double R = exp(r*(t/no_steps));
		double Rinv = 1.0/R;
		double u = exp(sigma*sqrt(t/no_steps));
		double d = 1.0/u;
		double uu= u*u;
		double pUp = (R-d)/(u-d);
		double pDown = 1.0 - pUp;
		vector<double> prices (no_steps+1);
		prices[0] = S*pow(d, no_steps);
		for (int i=1; i<=no_steps; ++i) prices[i] = uu*prices[i-1];
		vector<double> values (no_steps+1);
		for (int i=0; i<=no_steps; ++i) values[i] = generic_payoff(prices[i],K, PutCall, OpStyle);
		for (int CurrStep=no_steps-1 ; CurrStep>=1; --CurrStep) {
			for (int i=0; i<=CurrStep; ++i) {
				prices[i] = d*prices[i+1];
				values[i] = (pDown*values[i]+pUp*values[i+1])*Rinv;
			};
		}; 
		double delta = (values[1]-values[0])/(S*u-S*d);
		return delta;
	};


	static double CND( double X )
	{

		double L, K, w ;
		double const a1 = 0.31938153, a2 = -0.356563782, a3 = 1.781477937;
		double const a4 = -1.821255978, a5 = 1.330274429;

		L = fabs(X);
		K = 1.0 / (1.0 + 0.2316419 * L);
		w = 1.0 - 1.0 / sqrt(2 * Pi) * exp(-L *L / 2) * (a1 * K + a2 * K *K + a3 * pow(K,3) + a4 * pow(K,4) + a5 * pow(K,5));

		if (X < 0 ){
			w= 1.0 - w;
		}
		return w;
	}
	static double CountS(int n, double Spot, double K, double r, double v, double T, char PutCall, char OpStyle)
	{
		int i,j;
		vector<vector<double> >  S(n+1, vector<double>(n+1));
		vector<vector<double> > Op(n+1, vector<double>(n+1));
		double dt, u, d, p;
		dt = T/n;
		u = exp(v*sqrt(dt));
		d = 1/u;
		p = (exp(r*dt)-d) / (u-d);

		// Build the binomial tree
		for (j=0; j<=n; j++) {
			for (i=0; i<=j; i++) {
				S[i][j] = Spot*pow(u,j-i)*pow(d,i);
			}
		}

		// Compute terminal payoffs
		for (i=0; i<=n; i++) {
			if (PutCall=='C')
				Op[i][n] = max(S[i][n] - K, 0.0);
			else
				Op[i][n] = max(K - S[i][n], 0.0);
		}

		// Backward recursion through the tree
		for (j=n-1; j>=0; j--) {
			for (i=0; i<=j; i++) {
				if (OpStyle=='E')
					Op[i][j] = exp(-r*dt)*(p*(Op[i][j+1]) + (1-p)*(Op[i+1][j+1]));
				else {
					if (PutCall=='C')
						Op[i][j] = max(S[i][j] - K, exp(-r*dt)*(p*(Op[i][j+1]) + (1-p)*(Op[i+1][j+1])));
					else
						Op[i][j] = max(K - S[i][j], exp(-r*dt)*(p*(Op[i][j+1]) + (1-p)*(Op[i+1][j+1])));
				}
			}
		}
		return Op[0][0];
	}
	Binomial(void);
	~Binomial(void);
};

